In this script, we go through all the pre-registered proposed analyses. As a reminder, the main research questions where as follows:
We first present the main “sample theory based” analyses (also known as frequentist), separately on the North American and UK samples in parallel to answer our first two research questions, then on the total dataset to answer our third research question. We then provide additional Bayesian statistics where a null effect was found, as specified in the pre-registration.
# Library imports, general settings ==============
library(tidyverse); library(egg)
library(knitr)
library(lme4); library(lmerTest); library(simr)
# As in our discussion with Mike, we will use lmerTest for calculating p values
# in the mixed-effects models (noted as a deviation)
library(brms); library(rstan)
library(future)
plan(multicore, workers = 8) # Swap to multiprocess if running from Rstudio
library(lattice)
library(effects)
library(sjPlot)
library(robustlmm)
library(car)
# Load model comparison functions
source("helper/lrtests.R")
# Deal with package priority issues
select <- dplyr::select
theme_set(theme_bw(base_size = 10))
options("future" = T)
#knitr::opts_chunk$set(cache = TRUE)
print(sessionInfo()) #listing all info about R and packages info
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] furrr_0.2.1 future.apply_1.7.0 bridgesampling_1.0-0
## [4] car_3.0-10 robustlmm_2.4-4 sjPlot_2.8.7
## [7] effects_4.2-0 carData_3.0-4 lattice_0.20-41
## [10] future_1.21.0 rstan_2.21.2 StanHeaders_2.21.0-7
## [13] brms_2.14.4 Rcpp_1.0.5 simr_1.0.5
## [16] lmerTest_3.1-3 lme4_1.1-26 Matrix_1.3-2
## [19] knitr_1.30 egg_0.4.5 gridExtra_2.3
## [22] forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2
## [25] purrr_0.3.4 readr_1.4.0 tidyr_1.1.2
## [28] tibble_3.0.4 ggplot2_3.3.3 tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.1.0 htmlwidgets_1.5.3 grid_4.0.3
## [4] munsell_0.5.0 codetools_0.2-18 effectsize_0.4.1
## [7] statmod_1.4.35 DT_0.17 miniUI_0.1.1.1
## [10] withr_2.3.0 Brobdingnag_1.2-6 colorspace_2.0-0
## [13] rstudioapi_0.13 stats4_4.0.3 robustbase_0.93-7
## [16] bayesplot_1.8.0 listenv_0.8.0 emmeans_1.5.3
## [19] RLRsim_3.1-6 coda_0.19-4 parallelly_1.23.0
## [22] vctrs_0.3.6 generics_0.1.0 TH.data_1.0-10
## [25] xfun_0.20 R6_2.5.0 markdown_1.1
## [28] gamm4_0.2-6 projpred_2.0.2 assertthat_0.2.1
## [31] promises_1.1.1 scales_1.1.1 multcomp_1.4-15
## [34] nnet_7.3-14 gtable_0.3.0 globals_0.14.0
## [37] processx_3.4.5 sandwich_3.0-0 rlang_0.4.10
## [40] splines_4.0.3 broom_0.7.3 inline_0.3.17
## [43] yaml_2.2.1 reshape2_1.4.4 abind_1.4-5
## [46] modelr_0.1.8 threejs_0.3.3 crosstalk_1.1.1
## [49] backports_1.2.1 httpuv_1.5.5 rsconnect_0.8.16
## [52] tools_4.0.3 ellipsis_0.3.1 ggridges_0.5.3
## [55] plyr_1.8.6 base64enc_0.1-3 ps_1.5.0
## [58] prettyunits_1.1.1 zoo_1.8-8 haven_2.3.1
## [61] fs_1.5.0 survey_4.0 magrittr_2.0.1
## [64] data.table_1.13.6 openxlsx_4.2.3 colourpicker_1.1.0
## [67] reprex_0.3.0 mvtnorm_1.1-1 sjmisc_2.8.6
## [70] matrixStats_0.57.0 hms_1.0.0 shinyjs_2.0.0
## [73] mime_0.9 evaluate_0.14 xtable_1.8-4
## [76] shinystan_2.5.0 pbkrtest_0.5-0.1 rio_0.5.16
## [79] sjstats_0.18.1 readxl_1.3.1 ggeffects_1.0.1
## [82] rstantools_2.1.1 compiler_4.0.3 V8_3.4.0
## [85] crayon_1.3.4 minqa_1.2.4 htmltools_0.5.1
## [88] mgcv_1.8-33 later_1.1.0.1 RcppParallel_5.0.2
## [91] lubridate_1.7.9.2 DBI_1.1.0 sjlabelled_1.1.7
## [94] dbplyr_1.4.2 MASS_7.3-53 boot_1.3-25
## [97] cli_2.2.0 mitools_2.4 parallel_4.0.3
## [100] insight_0.11.1 igraph_1.2.6 pkgconfig_2.0.3
## [103] numDeriv_2016.8-1.1 foreign_0.8-81 binom_1.1-1
## [106] xml2_1.3.2 dygraphs_1.1.1.6 estimability_1.3
## [109] rvest_0.3.6 callr_3.5.1 digest_0.6.27
## [112] parameters_0.10.1 fastGHQuad_1.0 rmarkdown_2.6
## [115] cellranger_1.1.0 curl_4.3 shiny_1.5.0
## [118] gtools_3.8.2 nloptr_1.2.2.2 lifecycle_0.2.0
## [121] nlme_3.1-151 jsonlite_1.7.2 fansi_0.4.1
## [124] pillar_1.4.7 loo_2.4.1 fastmap_1.0.1
## [127] httr_1.4.2 plotrix_3.7-8 DEoptimR_1.0-8
## [130] pkgbuild_1.2.0 survival_3.2-7 glue_1.4.2
## [133] xts_0.12.1 bayestestR_0.8.0 zip_2.1.1
## [136] shinythemes_1.1.2 iterators_1.0.13 stringi_1.5.3
## [139] performance_0.6.1
# Read data ======================================
col_types <- cols(
labid = col_factor(),
subid = col_factor(),
subid_unique = col_factor(),
CDI.form = col_factor(),
CDI.nwords = col_integer(),
CDI.prop = col_number(),
CDI.agerange = col_factor(),
CDI.agedays = col_integer(),
CDI.agemin = col_integer(),
CDI.agemax = col_integer(),
vocab_nwords = col_integer(),
standardized.score.CDI = col_character(),
standardized.score.CDI.num = col_number(),
IDS_pref = col_number(),
language = col_factor(),
language_zone = col_factor(),
CDI.error = col_logical(),
Notes = col_character(),
trial_order = col_factor(),
method = col_factor(),
age_days = col_integer(),
age_mo = col_number(),
age_group = col_factor(),
nae = col_logical(),
gender = col_factor(),
second_session = col_logical()
)
data.total <- read_csv("data/02b_processed.csv", col_types = col_types)
# TODO: add saved results
Before moving on with the analysis, we have to ready the data by (a) checking for colinearity between z_age_months and CDI.z_age_months and correcting this if necessary, and (b) setting up the contrasts described in our data analysis.
First, we run a Kappa test on the possibility of colinearity between z_age_months and CDI.z_age_months.
# Run kappa test
k.age_months <- model.matrix(~ z_age_months + CDI.z_age_months, data = data.total) %>%
kappa(exact = T)
With a value of 3.1980728, we do not have a colinearity issue and can proceed with the data analysis as planned (The criteria of indicating colinearity is that kappa > 10).
We need gender as an effect-coded factor, and method as a deviation-coded factor. This is achieved in R by using the contr.sum() function with the number of levels for each factor. Notably, when subsetting the UK sample, only two levels of method out of the three in total were left.
# Set contrasts on the total dataset =============
contrasts(data.total$gender) <- contr.sum(2)
contrasts(data.total$method) <- contr.sum(3)
# Create sub-datasets, with contrasts ============
## NAE
data.nae <- data.total %>% subset(language_zone == "NAE") %>% droplevels()
contrasts(data.nae$gender) <- contr.sum(2)
contrasts(data.nae$method) <- contr.sum(3)
## UK
data.uk <- data.total %>% subset(language_zone == "British") %>% droplevels()
contrasts(data.uk$gender) <- contr.sum(2)
contrasts(data.uk$method) <- contr.sum(2) #note that UK sample has only 2 levels, so sum of zero contrasts set to 2 levels
## Other
data.other <- data.total %>% subset(language_zone == "Other") %>% droplevels()
contrasts(data.other$gender) <- contr.sum(2)
contrasts(data.other$method) <- contr.sum(3)
We first assess the amount of data we have overall per condition and their shape overall.
data.total %>%
group_by(language_zone, CDI.agerange, method, gender) %>%
summarise(N = n(), age = mean(CDI.agedays), sd = sd(CDI.agedays)) %>%
kable()
## `summarise()` regrouping output by 'language_zone', 'CDI.agerange', 'method' (override with `.groups` argument)
| language_zone | CDI.agerange | method | gender | N | age | sd |
|---|---|---|---|---|---|---|
| British | 18 | singlescreen | F | 24 | 551.2083 | 10.741950 |
| British | 18 | singlescreen | M | 24 | 550.5000 | 13.497182 |
| British | 18 | eyetracking | F | 9 | 548.5556 | 9.139353 |
| British | 18 | eyetracking | M | 11 | 547.8182 | 10.147100 |
| British | 24 | singlescreen | F | 18 | 741.6111 | 13.128948 |
| British | 24 | singlescreen | M | 15 | 745.0667 | 15.549307 |
| British | 24 | eyetracking | F | 13 | 731.5385 | 12.162890 |
| British | 24 | eyetracking | M | 5 | 737.8000 | 8.228001 |
| Other | 18 | singlescreen | F | 11 | 541.5455 | 7.160498 |
| Other | 18 | singlescreen | M | 13 | 544.3077 | 6.663140 |
| Other | 18 | eyetracking | F | 26 | 556.0769 | 14.133430 |
| Other | 18 | eyetracking | M | 30 | 560.8333 | 16.128970 |
| Other | 18 | hpp | F | 38 | 549.0526 | 10.812774 |
| Other | 18 | hpp | M | 38 | 555.1579 | 13.664961 |
| Other | 24 | singlescreen | F | 10 | 731.3000 | 14.802777 |
| Other | 24 | singlescreen | M | 12 | 721.1667 | 13.888081 |
| Other | 24 | eyetracking | F | 28 | 730.1786 | 9.996494 |
| Other | 24 | eyetracking | M | 25 | 730.1600 | 7.711896 |
| Other | 24 | hpp | F | 45 | 730.7333 | 17.357733 |
| Other | 24 | hpp | M | 35 | 730.5143 | 15.849237 |
| NAE | 18 | singlescreen | F | 19 | 554.9474 | 20.780530 |
| NAE | 18 | singlescreen | M | 14 | 581.2143 | 24.925052 |
| NAE | 18 | eyetracking | F | 1 | 549.0000 | NA |
| NAE | 18 | hpp | F | 64 | 557.9375 | 22.801333 |
| NAE | 18 | hpp | M | 66 | 556.1667 | 22.762035 |
| NAE | 24 | singlescreen | F | 23 | 737.1739 | 26.604258 |
| NAE | 24 | singlescreen | M | 20 | 739.6000 | 21.352678 |
| NAE | 24 | eyetracking | F | 2 | 756.5000 | 34.648232 |
| NAE | 24 | eyetracking | M | 1 | 745.0000 | NA |
| NAE | 24 | hpp | F | 58 | 731.6897 | 28.149480 |
| NAE | 24 | hpp | M | 65 | 726.2769 | 27.133068 |
More detailed information about Descriptive Statistics
#number of lab
data.total %>%
select(labid, language_zone) %>%
unique() %>%
group_by(language_zone) %>%
count()
## # A tibble: 3 x 2
## # Groups: language_zone [3]
## language_zone n
## <fct> <int>
## 1 British 4
## 2 Other 8
## 3 NAE 9
data.total %>%
group_by(language_zone, CDI.agerange) %>%
summarize(N = n())
## `summarise()` regrouping output by 'language_zone' (override with `.groups` argument)
## # A tibble: 6 x 3
## # Groups: language_zone [3]
## language_zone CDI.agerange N
## <fct> <fct> <int>
## 1 British 18 68
## 2 British 24 51
## 3 Other 18 156
## 4 Other 24 155
## 5 NAE 18 164
## 6 NAE 24 169
# age range in each age group and language zone
data.total %>%
select(subid, language_zone, CDI.agedays, CDI.agerange) %>%
unique() %>%
group_by(language_zone, CDI.agerange) %>%
summarize(age_min = (min(CDI.agedays)/365.25*12),
age_max = (max(CDI.agedays)/365.25*12))
## `summarise()` regrouping output by 'language_zone' (override with `.groups` argument)
## # A tibble: 6 x 4
## # Groups: language_zone [3]
## language_zone CDI.agerange age_min age_max
## <fct> <fct> <dbl> <dbl>
## 1 British 18 17.4 19.2
## 2 British 24 23.4 25.1
## 3 Other 18 17.4 19.5
## 4 Other 24 22.5 25.2
## 5 NAE 18 16.1 20.0
## 6 NAE 24 22.1 25.9
We then assess the data per lab in terms of sample size and CDI score (vocabulary size, for consistency between language zones).
by_lab <- data.total %>%
group_by(labid, language_zone, CDI.agerange) %>%
mutate(tested = n_distinct(subid_unique)) %>%
select(labid, language_zone, CDI.agerange, tested, vocab_nwords) %>%
nest(scores = vocab_nwords) %>%
mutate(model = map(scores, ~ lm(vocab_nwords ~ 1, data = .x)),
ci = map(model, confint)) %>%
transmute(tested = tested,
mean = map_dbl(model, ~ coefficients(.x)[[1]]),
ci_lower = map_dbl(ci, 1),
ci_upper = map_dbl(ci, 2)) %>%
arrange(language_zone) %>%
rownames_to_column()
# TODO: find a way to group by language zone?
ggplot(by_lab,
aes(x = labid, colour = language_zone,
y = mean, ymin = ci_lower, ymax = ci_upper)) +
geom_linerange() +
geom_point(aes(size = tested)) +
facet_grid(cols = vars(CDI.agerange), scales = "free") + coord_flip(ylim = c(0, 500)) +
xlab("Lab") + ylab("Vocabulary size") +
scale_colour_brewer(palette = "Dark2", name = "Language zone",
breaks = c("British", "NAE", "Other")) +
scale_size_continuous(name = "N", breaks = function(x) c(min(x), mean(x), max(x))) +
theme(legend.position = "bottom")
First, we want to assess quickly if there is a direct correlation between IDS preference and CDI score, computing a Pearson#’s product-moment correlation. We use standardized CDI scores for the North American sample, and raw scores for the British sample. [NOTE FROM MELANIE: I’m not sure this makes sense to do with the raw UK scores - since we’re collapsing across 18 and 24 month data and the data aren#’t standardized by age.].
# Statistics =====================================
## North American Sample
test.pearson.nae <- cor.test(data.nae$IDS_pref,
data.nae$z_standardized_CDI,
alternative = "two.sided", method = "pearson")
test.pearson.nae
##
## Pearson's product-moment correlation
##
## data: data.nae$IDS_pref and data.nae$z_standardized_CDI
## t = -0.91847, df = 305, p-value = 0.3591
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.16349833 0.05977293
## sample estimates:
## cor
## -0.05251901
## UK Sample
test.pearson.uk <- cor.test(data.uk$IDS_pref,
data.uk$z_vocab_nwords,
alternative = "two.sided", method = "pearson")
test.pearson.uk
##
## Pearson's product-moment correlation
##
## data: data.uk$IDS_pref and data.uk$z_vocab_nwords
## t = 1.0327, df = 117, p-value = 0.3039
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.08643398 0.27040987
## sample estimates:
## cor
## 0.09504018
Plots for correlation
## North American Sample
### Get correlation value for annotation
cor_text <- "paste(italic(R)^2, \" =\")"
cor_value <- round(test.pearson.nae$estimate, 3)
### Build plot
plot.pearson.nae <- data.nae %>%
ggplot(aes(x = IDS_pref,
y = standardized.score.CDI.num)) +
xlab("IDS preference") + ylab("Standardized CDI score") +
geom_point() +
geom_smooth(method = lm) +
annotate("text", x = -.9, y = 51, parse = T, size = 4,
label = paste(cor_text, cor_value, sep = "~"))
## UK Sample
cor_value <- round(test.pearson.uk$estimate, 3)
plot.pearson.uk <- data.uk %>%
ggplot(aes(x = IDS_pref,
y = vocab_nwords)) +
xlab("IDS preference") + ylab("Vocabulary size (in words)") +
geom_point() +
geom_smooth(method = lm) +
annotate("text", x = .8, y = 150, parse = T, size = 4,
label = paste(cor_text, cor_value, sep = "~"))
# Global plot
plot.pearson <- ggarrange(plot.pearson.nae, plot.pearson.uk, ncol = 2)
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 26 rows containing non-finite values (stat_smooth).
## Warning: Removed 26 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
plot.pearson
ggsave("plots/corr_plot.pdf", plot.pearson,
units = "mm", width = 180, height = 100, dpi = 1000)
We see no obvious direct link between IDS prefernce and CDI score here. However, an effect might appear once we take into account various factors that might interact with IDS preference and/or CDI score. We can also first enhance these plots with information about the age group at which infants were tested (18- or 24-month-old), using vocabulary size to better compare the NAE and UK samples.
plot.age_group <- data.total %>%
subset(language_zone != "Other") %>%
droplevels() %>%
ggplot(aes(x = IDS_pref,
y = vocab_nwords,
colour = CDI.agerange)) +
facet_wrap(vars(language_zone),
labeller = as_labeller(c("British" = "UK samples",
"NAE" = "North Amercian samples"))) +
xlab("Standardized IDS prefence score") + ylab("Vocabulary size (in words)") + theme(legend.position = "top") +
geom_point() +
geom_smooth(method = lm) +
scale_colour_brewer(palette = "Dark2", name = "Age group",
breaks = c("18", "24"),
labels = c("18mo", "24m"))
ggsave("plots/scatter_age.pdf", plot.age_group,
units = "mm", width = 180, height = 100, dpi = 1000)
## `geom_smooth()` using formula 'y ~ x'
(plot.age_group)
## `geom_smooth()` using formula 'y ~ x'
Here, we run a mixed-effects model including only theoretically motivated effects, as described in the pre-registration. We start with the full model bellow, simplifying the random effects structure until it converges.
# Run models =====================================
## NAE
data.nae$centered_IDS_pref <- scale(data.nae$IDS_pref, center = T, scale = F)
lmer.full.nae <- lmer(standardized.score.CDI.num ~ CDI.z_age_months + gender +
z_age_months + method + centered_IDS_pref +
centered_IDS_pref:method + centered_IDS_pref:CDI.z_age_months + centered_IDS_pref:z_age_months +
(1 | labid) + (1 | subid_unique),
data = data.nae)
summary(lmer.full.nae)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months +
## method + centered_IDS_pref + centered_IDS_pref:method + centered_IDS_pref:CDI.z_age_months +
## centered_IDS_pref:z_age_months + (1 | labid) + (1 | subid_unique)
## Data: data.nae
##
## REML criterion at convergence: 2806.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2240 -0.4550 -0.0619 0.4779 2.5238
##
## Random effects:
## Groups Name Variance Std.Dev.
## subid_unique (Intercept) 531.62 23.06
## labid (Intercept) 37.34 6.11
## Residual 246.81 15.71
## Number of obs: 307, groups: subid_unique, 211; labid, 8
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 32.47880 6.47774 49.79324 5.014
## CDI.z_age_months 0.04723 0.35956 121.44454 0.131
## gender1 -1.51119 1.88476 193.54899 -0.802
## z_age_months -0.01368 0.62999 69.41345 -0.022
## method1 6.62448 6.56943 122.84082 1.008
## method2 -13.97899 11.88500 177.54529 -1.176
## centered_IDS_pref -34.17431 24.81975 210.14326 -1.377
## method1:centered_IDS_pref 28.35879 25.46036 210.29032 1.114
## method2:centered_IDS_pref -59.58231 49.55101 209.34368 -1.202
## CDI.z_age_months:centered_IDS_pref 0.41746 1.10252 123.90420 0.379
## z_age_months:centered_IDS_pref 2.01470 2.07816 189.23346 0.969
## Pr(>|t|)
## (Intercept) 7.14e-06 ***
## CDI.z_age_months 0.896
## gender1 0.424
## z_age_months 0.983
## method1 0.315
## method2 0.241
## centered_IDS_pref 0.170
## method1:centered_IDS_pref 0.267
## method2:centered_IDS_pref 0.231
## CDI.z_age_months:centered_IDS_pref 0.706
## z_age_months:centered_IDS_pref 0.334
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) CDI.z__ gendr1 z_g_mn methd1 methd2 c_IDS_ m1:_ID m2:_ID
## CDI.z_g_mnt -0.058
## gender1 -0.031 0.022
## z_age_mnths 0.041 -0.037 -0.021
## method1 -0.716 0.000 0.021 0.105
## method2 0.880 -0.030 -0.035 0.017 -0.877
## cntrd_IDS_p 0.351 0.006 0.038 -0.040 -0.342 0.380
## mthd1:_IDS_ -0.330 -0.028 -0.027 0.017 0.356 -0.379 -0.927
## mthd2:_IDS_ 0.344 0.021 0.060 -0.008 -0.356 0.386 0.966 -0.971
## CDI.__:_IDS -0.011 0.012 -0.019 -0.019 -0.029 0.008 -0.049 0.027 -0.047
## z_g_m:_IDS_ -0.041 0.026 0.103 0.047 -0.037 0.022 0.028 -0.086 0.132
## CDI.__:
## CDI.z_g_mnt
## gender1
## z_age_mnths
## method1
## method2
## cntrd_IDS_p
## mthd1:_IDS_
## mthd2:_IDS_
## CDI.__:_IDS
## z_g_m:_IDS_ -0.071
# robust_lmer.full.nae <- robustlmm::rlmer(standardized.score.CDI.num ~ CDI.z_age_months + gender +
# z_age_months + method + centered_IDS_pref +
# centered_IDS_pref:method + centered_IDS_pref:CDI.z_age_months + centered_IDS_pref:z_age_months +
# (1 | labid),
# data = data.nae)
#
#
# summary(robust_lmer.full.nae) #this model is used to see if we can meet some statistical assumption better but we decided to use the original model as the inferential statistical results are consistent
full.nae_pvalue <- anova(lmer.full.nae) %>%
as_tibble(rownames = "Parameter") #this gives us the Type III p values
# ==========
# Sequentially removed random effects:
# IDS_pref:z_age_months
# IDS_pref:CDI.z_age_months
# IDS_pref:method
# IDS_pref
# method
# z_age_months
# gender
# CDI.z_age_months
#==========
#First, check linearity
# data.nae$resid <- residuals(lmer.full.nae)
#
# plot(data.nae$resid, data.nae$standardized.score.CDI)
#Second, check normality
plot_model(lmer.full.nae, type = 'diag') ## we do have right-skewed normality of residuals
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.18
## Current Matrix version is 1.3.2
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
## [[1]]
## `geom_smooth()` using formula 'y ~ x'
##
## [[2]]
## [[2]]$subid_unique
## `geom_smooth()` using formula 'y ~ x'
##
## [[2]]$labid
## `geom_smooth()` using formula 'y ~ x'
##
##
## [[3]]
##
## [[4]]
## `geom_smooth()` using formula 'y ~ x'
#Third, check autocorrelation
re_run_lme.full.nae <- nlme::lme(standardized.score.CDI.num ~ CDI.z_age_months + gender +
z_age_months + method + centered_IDS_pref +
centered_IDS_pref:method + centered_IDS_pref:CDI.z_age_months + centered_IDS_pref:z_age_months, random = ~1 | labid,
method = "REML",
data = data.nae, na.action = na.exclude)
plot(nlme::ACF(re_run_lme.full.nae, resType = "normalized")) #there is no sign for autocorrelation
#Lastly, check multi-collinearity
car::vif(lmer.full.nae) #we do see a multicollineartiy for the IDS preference variable, even though we have centered the IDS preference score. It is probably related to the number the participating labs (as this is the group level that we are controlling) and how we entered interaction between IDS preference and other variables (that lack variability in the current sample). We need to keep IDS preference in the model as exploring the relationship between IDS preference and CDI score is the key research question in the paper.
## GVIF Df GVIF^(1/(2*Df))
## CDI.z_age_months 1.009314 1 1.004646
## gender 1.039338 1 1.019479
## z_age_months 1.096811 1 1.047288
## method 1.284501 2 1.064593
## centered_IDS_pref 19.298596 1 4.393017
## method:centered_IDS_pref 20.886393 2 2.137794
## CDI.z_age_months:centered_IDS_pref 1.014832 1 1.007389
## z_age_months:centered_IDS_pref 1.292422 1 1.136847
lmer.full.uk <- lmer(vocab_nwords ~ CDI.z_age_months + gender +
z_age_months + method + IDS_pref +
IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
(1 | labid) + (1 | subid_unique),
data = data.uk)
summary(lmer.full.uk)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method +
## IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months +
## IDS_pref:z_age_months + (1 | labid) + (1 | subid_unique)
## Data: data.uk
##
## REML criterion at convergence: 1326
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.78175 -0.38449 0.01193 0.37856 1.94702
##
## Random effects:
## Groups Name Variance Std.Dev.
## subid_unique (Intercept) 5412.10 73.567
## labid (Intercept) 33.24 5.765
## Residual 2637.21 51.354
## Number of obs: 119, groups: subid_unique, 88; labid, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 174.753 11.544 3.754 15.138 0.000167 ***
## CDI.z_age_months 28.905 1.964 44.974 14.714 < 2e-16 ***
## gender1 19.043 9.820 78.517 1.939 0.056074 .
## z_age_months -5.065 3.455 70.444 -1.466 0.147086
## method1 -11.249 13.255 5.502 -0.849 0.431430
## IDS_pref 9.818 26.555 82.944 0.370 0.712523
## method1:IDS_pref -7.789 29.223 87.125 -0.267 0.790454
## CDI.z_age_months:IDS_pref 1.573 5.726 50.813 0.275 0.784675
## z_age_months:IDS_pref 1.154 8.805 89.735 0.131 0.896024
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) CDI.z__ gendr1 z_g_mn methd1 IDS_pr m1:IDS CDI.__:
## CDI.z_g_mnt 0.120
## gender1 0.044 -0.069
## z_age_mnths -0.011 -0.099 -0.019
## method1 -0.395 -0.074 -0.090 0.490
## IDS_pref -0.360 -0.074 -0.112 0.033 0.199
## mthd1:IDS_p 0.188 0.094 0.007 -0.117 -0.303 -0.203
## CDI.__:IDS_ -0.099 -0.328 0.056 0.064 0.094 0.090 -0.136
## z_g_mn:IDS_ -0.025 0.079 -0.244 -0.374 -0.099 -0.112 0.431 -0.228
full.uk_pvalue <- anova(lmer.full.uk) %>%
as_tibble(rownames = "Parameter") #this gives us the Type III p values
#==========
# Sequentially removed random effects:
# IDS_pref:z_age_months
# IDS_pref:CDI.z_age_months
# IDS_pref:method
# IDS_pref
# method
# z_age_months
# gender
# CDI.z_age_months
#First, check linearity. The plot looks linear
data.uk$resid <- residuals(lmer.full.uk)
plot(data.uk$resid, data.uk$vocab_nwords)
#Second, check normality
plot_model(lmer.full.uk, type = 'diag') ## we do have right-skewed normality of residuals
## [[1]]
## `geom_smooth()` using formula 'y ~ x'
##
## [[2]]
## [[2]]$subid_unique
## `geom_smooth()` using formula 'y ~ x'
##
## [[2]]$labid
## `geom_smooth()` using formula 'y ~ x'
##
##
## [[3]]
##
## [[4]]
## `geom_smooth()` using formula 'y ~ x'
#Third, check autocorrelation
re_run_lme.full.uk <- nlme::lme(vocab_nwords ~ CDI.z_age_months + gender +
z_age_months + method + IDS_pref +
IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months, random = ~1 | labid,
method = "REML",
data = data.nae, na.action = na.exclude)
plot(nlme::ACF(re_run_lme.full.uk, resType = "normalized")) #there is no sign for autocorrelation
#Lastly, check multi-collinearity
car::vif(lmer.full.uk) #no problem for multicollinearlity
## CDI.z_age_months gender z_age_months
## 1.142786 1.125370 1.678292
## method IDS_pref method:IDS_pref
## 1.592427 1.099496 1.457526
## CDI.z_age_months:IDS_pref z_age_months:IDS_pref
## 1.188763 1.737367
We now want to check the statistical power of significant effects, and discard any models with significant effects that do not reach 80% power. This however leads to too many warnings of singularity issues on the model updates inherent to the simr power simulations, hence we cannot obtain satisfactory power estimates as pre-registered.
AST: Note that we don’t have any IV(s) that turned out to be significant in the Full NAE model. So we won’t run the power analysis check. For the UK full model, there are two statistically significant IV: CDI_age and gender. The post hoc power check suggested that we have high power in detecting the effect of CDI_age but not gender. Note that gender has a smaller effect size to begin with, so this may partially explain why we have less power in detecting it in the model. As there can be a number of different factors that determines the posthoc power, we decided not to remove gender in the model based on posthoc power analysis check.
check_pwr_uk_cdi_age <- simr::powerSim(lmer.full.uk, test = fixed("CDI.z_age_months", method = "z"), seed = 2, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into
check_pwr_uk_cdi_age
check_pwr_uk_gender <- simr::powerSim(lmer.full.uk, test = fixed("gender", method = "z"), seed = 3, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into
check_pwr_uk_gender
For this combined analysis, we first need to restrain the age range for the NAE sample (previously ±2 months, now ±1.5 months).
# Create dataset with British and NAE only
data.uk_nae <- data.total %>%
subset(language_zone %in% c("British", "NAE")) %>%
mutate(CDI.agemin = ifelse(language_zone == "NAE",
CDI.agemin + round(.5*365.25/12),
CDI.agemin),
CDI.agemax = ifelse(language_zone == "NAE",
CDI.agemax - round(.5*365.25/12),
CDI.agemax)) %>%
subset(!(CDI.agedays < CDI.agemin | CDI.agedays > CDI.agemax)) %>%
droplevels()
# Create contrasts for analysis
contrasts(data.uk_nae$gender) <- contr.sum(2)
contrasts(data.uk_nae$method) <- contr.sum(3)
contrasts(data.uk_nae$language_zone) <- contr.sum(2)
We can then run the planned combined analysis adding the main effect and interactions of language_zone.
lmer.full.uk_nae <- lmer(CDI.prop ~ CDI.z_age_months + language_zone + gender +
z_age_months + method + IDS_pref + IDS_pref:language_zone +
IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
(1 + CDI.z_age_months | labid) + (1 | subid_unique),
data = data.uk_nae)
summary(lmer.full.uk_nae)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CDI.prop ~ CDI.z_age_months + language_zone + gender + z_age_months +
## method + IDS_pref + IDS_pref:language_zone + IDS_pref:method +
## IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 +
## CDI.z_age_months | labid) + (1 | subid_unique)
## Data: data.uk_nae
##
## REML criterion at convergence: -61.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.83614 -0.46389 -0.04311 0.41335 2.60155
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subid_unique (Intercept) 0.0267932 0.16369
## labid (Intercept) 0.0003603 0.01898
## CDI.z_age_months 0.0009806 0.03131 0.44
## Residual 0.0184535 0.13584
## Number of obs: 424, groups: subid_unique, 302; labid, 13
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.320830 0.017878 12.110147 17.945 4.31e-10 ***
## CDI.z_age_months 0.043700 0.009272 13.321953 4.713 0.00038 ***
## language_zone1 0.087738 0.021795 10.500676 4.026 0.00219 **
## gender1 0.032723 0.011939 270.491415 2.741 0.00654 **
## z_age_months -0.002774 0.004295 164.183977 -0.646 0.51919
## method1 -0.007741 0.024041 21.324230 -0.322 0.75061
## method2 0.009802 0.035669 15.572064 0.275 0.78708
## IDS_pref 0.032520 0.041549 295.409651 0.783 0.43445
## language_zone1:IDS_pref 0.036505 0.056770 318.273301 0.643 0.52066
## method1:IDS_pref -0.040507 0.053921 317.820649 -0.751 0.45307
## method2:IDS_pref -0.018199 0.087818 306.828123 -0.207 0.83596
## CDI.z_age_months:IDS_pref 0.009728 0.008057 172.834552 1.207 0.22889
## z_age_months:IDS_pref -0.001722 0.011708 278.006432 -0.147 0.88318
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
combined.full.uk_nae_pvalue <- anova(lmer.full.uk_nae) %>%
as_tibble(rownames = "Parameter") #this gives us the Type III p values
#==========
# Sequentially removed random effects:
# IDS_pref:z_age_months
# IDS_pref:CDI.z_age_months
# IDS_pref:method
# IDS_pref:language_zone
# IDS_pref
# method
# z_age_months
# gender
# language_zone
#==========
#First, check linearity. The plot looks linear
data.uk_nae$resid <- residuals(lmer.full.uk_nae)
plot(data.uk_nae$resid, data.uk_nae$CDI.prop)
#Second, check normality
plot_model(lmer.full.uk_nae, type = 'diag') ## we do have right-skewed normality of residuals
## [[1]]
## `geom_smooth()` using formula 'y ~ x'
##
## [[2]]
## [[2]]$subid_unique
## `geom_smooth()` using formula 'y ~ x'
##
## [[2]]$labid
## `geom_smooth()` using formula 'y ~ x'
##
##
## [[3]]
##
## [[4]]
## `geom_smooth()` using formula 'y ~ x'
#Third, check autocorrelation
re_run_lme.full.uk_nae <- nlme::lme(CDI.prop ~ CDI.z_age_months + language_zone + gender +
z_age_months + method + IDS_pref + IDS_pref:language_zone +
IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months,
random = ~CDI.z_age_months | labid,
method = "REML",
data = data.uk_nae, na.action = na.exclude)
plot(nlme::ACF(re_run_lme.full.uk_nae, resType = "normalized")) #there is no sign for autocorrelation
#Lastly, check multi-collinearity
car::vif(lmer.full.uk_nae) #no problem for multicollinearlity
## GVIF Df GVIF^(1/(2*Df))
## CDI.z_age_months 1.019950 1 1.009926
## language_zone 2.325115 1 1.524833
## gender 1.036435 1 1.018055
## z_age_months 1.474609 1 1.214335
## method 3.284643 2 1.346239
## IDS_pref 1.447535 1 1.203135
## language_zone:IDS_pref 2.971184 1 1.723712
## method:IDS_pref 3.985056 2 1.412891
## CDI.z_age_months:IDS_pref 1.053846 1 1.026570
## z_age_months:IDS_pref 1.430296 1 1.195950
We then compute \(p\)-values, but leave out power estimates for those \(p\)-values as above. Again, we have a lot of singular fit issues for the power checks and decided not to remove parameters based on posthoc power analysis.
check_pwr_combined_cdi_age <- simr::powerSim(lmer.full.uk_nae, test = fixed("CDI.z_age_months", method = "z"), seed = 2, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into
check_pwr_combined_cdi_age
check_pwr_combined_lang_zone <- simr::powerSim(lmer.full.uk_nae, test = fixed("language_zone", method = "z"), seed = 3, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into
check_pwr_combined_lang_zone
check_pwr_combined_gender <- simr::powerSim(lmer.full.uk_nae, test = fixed("gender", method = "z"), seed = 3, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into
check_pwr_combined_gender
Given that the `frequentist#’ tests above did not discover any expected significant effects, we now run a Bayesian analysis, specifically Bayes factors (or \(b\)-values) from model comparisons, to determine whether our data provide evidence for a true null effect or merely fail to provide evidence either way.
We first run models on the separate UK and NAE samples.
First of all, we need to set up the priors for the Bayeisan analysis. In the RR, we said that we would use a Cauchy distribution for priors over fixed effects with normally distributed priors over random effect parameters."
Set up a really weak and not too informative prior
prior <-c(set_prior("normal(10, 1)", class = "Intercept"), #we can't have negative CDI
set_prior("cauchy(0, 1)", class = "b"), #all IVs have the same prior..
set_prior("normal(0, 1)", class = "sd"))
prior_combined_mod <-c(set_prior("normal(10, 1)", class = "Intercept"),
set_prior("cauchy(0, 1)", class = "b"), #all IVs have the same prior..
set_prior("normal(0, 1)", class = "sd"),
set_prior("lkj(2)", class = "L")) #adding this as there is a random slope and we need a correlation matrix prior between slope and intercept
Set up another prior that mainly changed the meaan for the fixed effect and see if there is big difference? Note: numbers changed but not to the extent that change BF conclusion.
prior_big <-c(set_prior("normal(10, 1)", class = "Intercept"), #we can't have negative CDI
set_prior("cauchy(10, 5)", class = "b"), #all IVs have the same prior..
set_prior("normal(0, 1)", class = "sd"))
prior_combined_big <-c(set_prior("normal(10, 1)", class = "Intercept"),
set_prior("cauchy(10, 5)", class = "b"), #all IVs have the same prior..
set_prior("normal(0, 1)", class = "sd"),
set_prior("lkj(2)", class = "L")) #adding this as there is a random slope and we need a correlation matrix prior between slope and intercept
We also try more informative prior.
prior_nae_1 <-c(set_prior("normal(50, 10)", class = "Intercept"), #note that nae's DV is percentile score, but I will assume at 18 months, the average is at 50 percentile
set_prior("cauchy(5, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 5 % percentile increase?
set_prior("cauchy(-15, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males
set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
set_prior("cauchy(0, 1)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
set_prior("normal(0, 1)", class = "sd"))
prior_nae_null_1 <-c(set_prior("normal(50, 10)", class = "Intercept"), #note that nae's DV is percentile score, but I will assume at 18 months, the average is at 50 percentile
set_prior("cauchy(5, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 5 % percentile increase?
set_prior("cauchy(-15, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males
set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
set_prior("cauchy(0, 1)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
set_prior("normal(0, 1)", class = "sd"))
prior_nae_2 <-c(set_prior("normal(50, 10)", class = "Intercept"), #note that nae's DV is percentile score, but I will assume at 18 months, the average is at 50 percentile
set_prior("cauchy(5, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 5 % percentile increase?
set_prior("cauchy(-15, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males
set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
set_prior("cauchy(0, 1)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
#note all interaction priors are weak as we don't really have much info to form a prior
set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
set_prior("cauchy(0, 1)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
set_prior("cauchy(0, 1)", class = "b", coef = "method1:IDS_pref"),
set_prior("cauchy(0, 1)", class = "b", coef = "method2:IDS_pref"),
set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months:IDS_pref"),
set_prior("normal(0, 1)", class = "sd"))
prior_nae_null_2 <-c(set_prior("normal(50, 10)", class = "Intercept"), #note that nae's DV is percentile score, but I will assume at 18 months, the average is at 50 percentile
set_prior("cauchy(5, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 5 % percentile increase?
set_prior("cauchy(-15, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males
set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
set_prior("cauchy(0, 1)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
#note all interaction priors are weak as we don't really have much info to form a prior
set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
set_prior("cauchy(0, 1)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
set_prior("cauchy(0, 1)", class = "b", coef = "method1:IDS_pref"),
set_prior("cauchy(0, 1)", class = "b", coef = "method2:IDS_pref"),
set_prior("normal(0, 1)", class = "sd"))
prior_uk_1 <-c(set_prior("normal(40, 20)", class = "Intercept"), #note that UK's DV is raw score, based on wordbank (oxford CDI), the average (50 percentile) is around 41.6
set_prior("cauchy(15, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase (proxy from wordbank)
set_prior("cauchy(-10, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males
set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
set_prior("normal(0, 1)", class = "sd"))
prior_uk_null_1 <-c(set_prior("normal(40, 20)", class = "Intercept"), #note that UK's DV is raw score, based on wordbank (oxford CDI), the average (50 percentile) is around 41.6
set_prior("cauchy(15, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase (proxy from wordbank)
set_prior("cauchy(-10, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males
set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
set_prior("normal(0, 1)", class = "sd"))
prior_uk_2 <-c(set_prior("normal(40, 20)", class = "Intercept"), #note that UK's DV is raw score, based on wordbank (oxford CDI), the average (50 percentile) is around 41.6
set_prior("cauchy(15, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase (proxy from wordbank)
set_prior("cauchy(-10, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males
set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
#note all interaction priors are weak as we don't really have much info to form a prior
set_prior("cauchy(0, 1)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
set_prior("cauchy(0, 1)", class = "b", coef = "method1:IDS_pref"),
set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months:IDS_pref"),
set_prior("normal(0, 1)", class = "sd"))
prior_uk_null_2 <-c(set_prior("normal(40, 20)", class = "Intercept"), #note that UK's DV is raw score, based on wordbank (oxford CDI), the average (50 percentile) is around 41.6
set_prior("cauchy(15, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase (proxy from wordbank)
set_prior("cauchy(-10, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males
set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
#note all interaction priors are weak as we don't really have much info to form a prior
set_prior("cauchy(0, 1)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
set_prior("cauchy(0, 1)", class = "b", coef = "method1:IDS_pref"),
set_prior("normal(0, 1)", class = "sd"))
# Note that NAE CDI WS total score is 680 and UK CDI total score is 416.
# NAE average scores for 50th percentile is 76.3, so it is 76.3/680 ~ 11% and UK CDI average scores for 50th percentile is 41.6, so it is 41.6/416 ~ 10%
prior_combined_full_info <-c(set_prior("normal(0.1, 0.05)", class = "Intercept"), # given the above info, we expect on average, proportion score is 10%
set_prior("cauchy(0.03, 0.01)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase in UK cdi and lead to 28 words increase in NAE (proxy from wordbank). In general, 15/416 = 3% and 28/680 = 4%
set_prior("cauchy(-0.03, 0.01)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males, 11/416 ~ 2.6% and Wordbank NAE suggested females knows approximately 27 words than males, 27/680 ~ 3.9%
set_prior("cauchy(0, 0.01)", class = "b", coef = "language_zone1"), #weak prior as don't know the effect of language zone to CDI score
set_prior("cauchy(0, 0.01)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
set_prior("cauchy(0, 0.01)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
set_prior("cauchy(0, 0.01)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
set_prior("cauchy(0, 0.01)", class = "b", coef = "z_age_months"),
#note all interaction priors are weak as we don't really have much info to form a prior
set_prior("cauchy(0, 0.01)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
set_prior("cauchy(0, 0.01)", class = "b", coef = "language_zone1:IDS_pref"),
set_prior("cauchy(0, 0.01)", class = "b", coef = "method1:IDS_pref"),
set_prior("cauchy(0, 0.01)", class = "b", coef = "method2:IDS_pref"),
set_prior("cauchy(0, 0.01)", class = "b", coef = "z_age_months:IDS_pref"),
set_prior("normal(0, 0.01)", class = "sd"),
set_prior("lkj(2)", class = "L")) #adding this as there is a random slope and we need a correlation matrix prior between slope and intercept
prior_combined_null_info <-c(set_prior("normal(0.1, 0.05)", class = "Intercept"), # given the above info, we expect on average, proportion score is 10%
set_prior("cauchy(0.03, 0.01)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase in UK cdi and lead to 28 words increase in NAE (proxy from wordbank). In general, 15/416 = 3% and 28/680 = 4%
set_prior("cauchy(-0.03, 0.01)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males, 11/416 ~ 2.6% and Wordbank NAE suggested females knows approximately 27 words than males, 27/680 ~ 3.9%
set_prior("cauchy(0, 0.01)", class = "b", coef = "language_zone1"), #weak prior as don't know the effect of language zone to CDI score
set_prior("cauchy(0, 0.01)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
set_prior("cauchy(0, 0.01)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
set_prior("cauchy(0, 0.01)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
set_prior("cauchy(0, 0.01)", class = "b", coef = "z_age_months"),
#note all interaction priors are weak as we don't really have much info to form a prior
set_prior("cauchy(0, 0.01)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
set_prior("cauchy(0, 0.01)", class = "b", coef = "method1:IDS_pref"),
set_prior("cauchy(0, 0.01)", class = "b", coef = "method2:IDS_pref"),
set_prior("cauchy(0, 0.01)", class = "b", coef = "z_age_months:IDS_pref"),
set_prior("normal(0, 0.01)", class = "sd"),
set_prior("lkj(2)", class = "L")) #adding this as there is a random slope and we need a correlation matrix prior between slope and intercept
In the following analysis, we follow what we have proposed in the RR about calculating the Bayes factor for the three variables of interest: 1) The main effect IDS preference on CDI scores 2) The interaction effect between IDS preference and IDS test age on CDI scores 3) The interaction effect between language_zone (i.e., dialect in RR) and IDS preference on CDI scores
For the first two variables, we run it on the NAE and UK model separately. For the third variable, we run in on the combined UK & NAE model.
To calcualte Bayes factor for the main effect IDS preference on CDI scores, we need to two models: a full model that contains the main effect of variables of interest, and the other null model that does not contain the main effect of interest. As all the interaction terms in the models contains IDS preference, in the following, we only include main effects in the full and null models for comparisons.
Sys.setenv(R_FUTURE_RNG_ONMISUSE = "ignore") #this line of code is needed to get rid of warning message due to the "future" package. https://github.com/satijalab/seurat/issues/3622
lmer.bf_full_1.nae <- brm(standardized.score.CDI.num ~ CDI.z_age_months + gender +
z_age_months + method + IDS_pref +
(1 | labid) + (1 | subid_unique),
data = data.nae,
family = gaussian,
prior = prior_nae_1,
iter = 4000,
warmup = 2000,
chains = 4,
save_pars = save_pars(all = T),
control = list(adapt_delta = .95),
seed=456)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
##
## SAMPLING FOR MODEL '6361f04d0b2eda23c05b2e4be9a58962' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 7.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.75 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 2.19751 seconds (Warm-up)
## Chain 1: 2.05293 seconds (Sampling)
## Chain 1: 4.25045 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '6361f04d0b2eda23c05b2e4be9a58962' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 5.9e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.59 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 2.35878 seconds (Warm-up)
## Chain 2: 2.07516 seconds (Sampling)
## Chain 2: 4.43394 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '6361f04d0b2eda23c05b2e4be9a58962' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 6.4e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.64 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 2.25339 seconds (Warm-up)
## Chain 3: 1.3515 seconds (Sampling)
## Chain 3: 3.60489 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '6361f04d0b2eda23c05b2e4be9a58962' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 6.2e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.62 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 2.31966 seconds (Warm-up)
## Chain 4: 2.0556 seconds (Sampling)
## Chain 4: 4.37526 seconds (Total)
## Chain 4:
summary(lmer.bf_full_1.nae)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + (1 | labid) + (1 | subid_unique)
## Data: data.nae (Number of observations: 307)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~labid (Number of levels: 8)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.05 0.75 0.05 2.78 1.00 5477 4927
##
## ~subid_unique (Number of levels: 211)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.87 0.66 0.03 2.44 1.00 6047 4451
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 40.24 2.10 36.10 44.27 1.00 11676 5899
## CDI.z_age_months 0.14 0.55 -0.92 1.20 1.00 16386 5805
## gender1 -1.98 1.61 -5.11 1.19 1.00 15997 5403
## z_age_months -0.09 0.43 -0.93 0.75 1.00 14152 5721
## method1 0.19 1.28 -2.31 2.96 1.00 10776 4501
## method2 -0.05 1.79 -3.86 3.55 1.00 8843 3907
## IDS_pref -0.62 2.16 -6.38 2.91 1.00 7741 3066
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 28.08 1.17 25.89 30.48 1.00 13999 5427
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_1.nae <- brm(standardized.score.CDI.num ~ CDI.z_age_months + gender +
z_age_months + method +
(1 | labid) + (1 | subid_unique),
data = data.nae,
family = gaussian,
prior = prior_nae_null_1,
iter = 4000,
warmup = 2000,
chains = 4,
save_pars = save_pars(all = T),
control = list(adapt_delta = .95),
seed =123)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
##
## SAMPLING FOR MODEL '4bb864c3fbd97441232e09115f930fe8' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 8.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.81 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 2.40827 seconds (Warm-up)
## Chain 1: 1.76915 seconds (Sampling)
## Chain 1: 4.17742 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '4bb864c3fbd97441232e09115f930fe8' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.00012 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.2 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 2.0687 seconds (Warm-up)
## Chain 2: 2.04309 seconds (Sampling)
## Chain 2: 4.11179 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '4bb864c3fbd97441232e09115f930fe8' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.000298 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 2.98 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 2.14399 seconds (Warm-up)
## Chain 3: 2.02593 seconds (Sampling)
## Chain 3: 4.16992 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '4bb864c3fbd97441232e09115f930fe8' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 7.9e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.79 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 2.128 seconds (Warm-up)
## Chain 4: 2.0381 seconds (Sampling)
## Chain 4: 4.1661 seconds (Total)
## Chain 4:
summary(lmer.bf_null_1.nae)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months + method + (1 | labid) + (1 | subid_unique)
## Data: data.nae (Number of observations: 307)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~labid (Number of levels: 8)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.05 0.74 0.05 2.74 1.00 6284 4654
##
## ~subid_unique (Number of levels: 211)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.86 0.66 0.03 2.42 1.00 5865 4466
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 40.14 2.11 35.93 44.29 1.00 13214 5661
## CDI.z_age_months 0.14 0.55 -0.93 1.22 1.00 16648 5387
## gender1 -1.97 1.65 -5.17 1.24 1.00 17327 5420
## z_age_months -0.11 0.42 -0.95 0.71 1.00 16529 5611
## method1 0.22 1.22 -2.15 2.89 1.00 10974 4928
## method2 -0.09 1.81 -4.07 3.51 1.00 9396 3756
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 28.10 1.15 25.97 30.42 1.00 17149 5479
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Calculate the marginal log likelihoods for each hypothese using Bridge sampling
marloglik_full_1.nae <- bridge_sampler(lmer.bf_full_1.nae, silent = T)
marloglik_null_1.nae <- bridge_sampler(lmer.bf_null_1.nae, silent = T)
Bayes factor
BF_full_1.nae <- bayes_factor(marloglik_full_1.nae, marloglik_null_1.nae)
BF_full_1.nae
## Estimated Bayes factor in favor of x1 over x2: 0.63608
lmer.bf_full_2.nae <- brm(standardized.score.CDI.num ~ CDI.z_age_months + gender +
z_age_months + method + IDS_pref +
IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
(1 | labid) + (1 | subid_unique),
data = data.nae,
family = gaussian,
prior = prior_nae_2,
iter = 4000,
warmup = 2000,
chains = 4,
save_pars = save_pars(all = T),
control = list(adapt_delta = .95),
seed=890)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
##
## SAMPLING FOR MODEL 'f3788e4416fd7ee1c5e7b74b10adc8b6' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000157 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.57 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 2.48097 seconds (Warm-up)
## Chain 1: 2.15046 seconds (Sampling)
## Chain 1: 4.63143 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'f3788e4416fd7ee1c5e7b74b10adc8b6' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.000325 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 3.25 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 2.38087 seconds (Warm-up)
## Chain 2: 2.13843 seconds (Sampling)
## Chain 2: 4.5193 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'f3788e4416fd7ee1c5e7b74b10adc8b6' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.000124 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.24 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 2.27497 seconds (Warm-up)
## Chain 3: 2.16732 seconds (Sampling)
## Chain 3: 4.44229 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'f3788e4416fd7ee1c5e7b74b10adc8b6' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 8.8e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.88 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 2.49946 seconds (Warm-up)
## Chain 4: 2.14814 seconds (Sampling)
## Chain 4: 4.64759 seconds (Total)
## Chain 4:
summary(lmer.bf_full_2.nae)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 | labid) + (1 | subid_unique)
## Data: data.nae (Number of observations: 307)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~labid (Number of levels: 8)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.03 0.75 0.04 2.78 1.00 5677 3586
##
## ~subid_unique (Number of levels: 211)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.89 0.66 0.04 2.46 1.00 7136 4836
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 40.08 2.23 35.69 44.40 1.00 8169
## CDI.z_age_months 0.08 0.58 -1.09 1.23 1.00 18071
## gender1 -1.91 1.61 -5.03 1.27 1.00 19483
## z_age_months -0.18 0.44 -1.07 0.68 1.00 18190
## method1 0.23 1.38 -2.36 3.19 1.00 7189
## method2 -0.16 2.16 -4.71 3.72 1.00 5468
## IDS_pref -0.91 2.51 -7.88 2.79 1.00 7484
## method1:IDS_pref -0.01 2.04 -4.54 4.43 1.00 11545
## method2:IDS_pref 0.35 2.54 -4.30 6.75 1.00 7882
## CDI.z_age_months:IDS_pref 0.45 1.07 -1.44 2.90 1.00 13628
## z_age_months:IDS_pref 0.81 1.15 -1.11 3.49 1.00 13396
## Tail_ESS
## Intercept 4067
## CDI.z_age_months 5831
## gender1 5281
## z_age_months 5450
## method1 2950
## method2 2415
## IDS_pref 3357
## method1:IDS_pref 3582
## method2:IDS_pref 2938
## CDI.z_age_months:IDS_pref 5707
## z_age_months:IDS_pref 5809
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 28.09 1.17 25.92 30.45 1.00 19893 5820
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_2.nae <- brm(standardized.score.CDI.num ~ CDI.z_age_months + gender +
z_age_months + method + IDS_pref +
IDS_pref:method + IDS_pref:CDI.z_age_months +
(1 | labid) + (1 | subid_unique),
data = data.nae,
family = gaussian,
prior = prior_nae_null_2,
iter = 4000,
warmup = 2000,
chains = 4,
save_pars = save_pars(all = T),
control = list(adapt_delta = .95),
seed =102)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
##
## SAMPLING FOR MODEL '0e2a1080bab56c0e90a24e05388bf5a7' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000255 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.55 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 2.42628 seconds (Warm-up)
## Chain 1: 2.09945 seconds (Sampling)
## Chain 1: 4.52572 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '0e2a1080bab56c0e90a24e05388bf5a7' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.000108 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.08 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 2.28607 seconds (Warm-up)
## Chain 2: 2.11654 seconds (Sampling)
## Chain 2: 4.40261 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '0e2a1080bab56c0e90a24e05388bf5a7' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.000263 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 2.63 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 2.48995 seconds (Warm-up)
## Chain 3: 2.13717 seconds (Sampling)
## Chain 3: 4.62712 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '0e2a1080bab56c0e90a24e05388bf5a7' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.000267 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 2.67 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 2.32855 seconds (Warm-up)
## Chain 4: 2.10664 seconds (Sampling)
## Chain 4: 4.43519 seconds (Total)
## Chain 4:
summary(lmer.bf_null_2.nae)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + (1 | labid) + (1 | subid_unique)
## Data: data.nae (Number of observations: 307)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~labid (Number of levels: 8)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.04 0.75 0.04 2.71 1.00 5908 4221
##
## ~subid_unique (Number of levels: 211)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.88 0.65 0.04 2.41 1.00 6257 4689
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 40.17 2.10 35.97 44.34 1.00 12196
## CDI.z_age_months 0.08 0.57 -1.04 1.21 1.00 17288
## gender1 -1.97 1.60 -5.09 1.16 1.00 20629
## z_age_months -0.10 0.42 -0.95 0.75 1.00 17914
## method1 0.21 1.29 -2.22 3.01 1.00 10045
## method2 -0.10 1.90 -4.22 3.83 1.00 8727
## IDS_pref -0.61 2.16 -6.38 2.75 1.00 8840
## method1:IDS_pref -0.07 2.00 -4.55 4.17 1.00 12060
## method2:IDS_pref 0.14 2.36 -4.77 5.66 1.00 8495
## CDI.z_age_months:IDS_pref 0.47 1.06 -1.44 2.86 1.00 15876
## Tail_ESS
## Intercept 5688
## CDI.z_age_months 6153
## gender1 6037
## z_age_months 5889
## method1 4185
## method2 3640
## IDS_pref 3328
## method1:IDS_pref 3278
## method2:IDS_pref 3061
## CDI.z_age_months:IDS_pref 5247
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 28.11 1.15 25.96 30.42 1.00 19020 6352
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Calculate the marginal log likelihoods for each hypothese using Bridge sampling
marloglik_full_2.nae <- bridge_sampler(lmer.bf_full_2.nae, silent = T)
marloglik_null_2.nae <- bridge_sampler(lmer.bf_null_2.nae, silent = T)
Bayes factor
BF_full_2.nae <- bayes_factor(marloglik_full_2.nae, marloglik_null_2.nae)
BF_full_2.nae
## Estimated Bayes factor in favor of x1 over x2: 0.98500
lmer.bf_full_1.uk <- brm(vocab_nwords ~ CDI.z_age_months + gender +
z_age_months + method + IDS_pref +
(1 | labid) + (1 | subid_unique),
data = data.uk,
family = gaussian,
prior = prior_uk_1,
iter = 4000,
warmup = 2000,
chains = 4,
save_pars = save_pars(all = T),
control = list(adapt_delta = .95),
seed=213)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
##
## SAMPLING FOR MODEL '379ff944c7b93801a64060c542798150' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 7e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.7 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 1.17577 seconds (Warm-up)
## Chain 1: 0.925538 seconds (Sampling)
## Chain 1: 2.10131 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '379ff944c7b93801a64060c542798150' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7.9e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.79 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.896499 seconds (Warm-up)
## Chain 2: 0.956413 seconds (Sampling)
## Chain 2: 1.85291 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '379ff944c7b93801a64060c542798150' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 7e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.7 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 1.18867 seconds (Warm-up)
## Chain 3: 0.946536 seconds (Sampling)
## Chain 3: 2.1352 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '379ff944c7b93801a64060c542798150' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 6.1e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.61 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 1.05547 seconds (Warm-up)
## Chain 4: 1.53038 seconds (Sampling)
## Chain 4: 2.58585 seconds (Total)
## Chain 4:
summary(lmer.bf_full_1.uk)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + (1 | labid) + (1 | subid_unique)
## Data: data.uk (Number of observations: 119)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~labid (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.81 0.61 0.04 2.28 1.00 5852 3538
##
## ~subid_unique (Number of levels: 88)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.79 0.61 0.03 2.27 1.00 5496 3000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 153.94 8.21 137.27 169.56 1.00 9730 5591
## CDI.z_age_months 27.35 2.94 21.38 32.85 1.00 12045 5528
## gender1 9.90 9.69 -8.52 28.64 1.00 10527 5291
## z_age_months -1.03 1.63 -5.02 1.55 1.00 8414 4845
## method1 -0.69 3.17 -9.13 4.12 1.00 5119 1963
## IDS_pref 0.59 5.30 -7.24 13.01 1.00 4256 1457
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 91.43 6.38 80.20 105.07 1.00 9396 5966
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_1.uk <- brm(vocab_nwords ~ CDI.z_age_months + gender +
z_age_months + method +
(1 | labid) + (1 | subid_unique),
data = data.uk,
family = gaussian,
prior = prior_uk_null_1,
iter = 4000,
warmup = 2000,
chains = 4,
save_pars = save_pars(all = T),
control = list(adapt_delta = .95),
seed =312)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
##
## SAMPLING FOR MODEL '844f09eb03d0f7d774bc9ac00a910248' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.00023 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.3 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.917603 seconds (Warm-up)
## Chain 1: 0.853581 seconds (Sampling)
## Chain 1: 1.77118 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '844f09eb03d0f7d774bc9ac00a910248' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 4.1e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.41 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.769098 seconds (Warm-up)
## Chain 2: 0.95405 seconds (Sampling)
## Chain 2: 1.72315 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '844f09eb03d0f7d774bc9ac00a910248' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 6.2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.62 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.75625 seconds (Warm-up)
## Chain 3: 0.94152 seconds (Sampling)
## Chain 3: 1.69777 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '844f09eb03d0f7d774bc9ac00a910248' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.000204 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 2.04 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.77893 seconds (Warm-up)
## Chain 4: 0.942255 seconds (Sampling)
## Chain 4: 1.72118 seconds (Total)
## Chain 4:
summary(lmer.bf_null_1.uk)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method + (1 | labid) + (1 | subid_unique)
## Data: data.uk (Number of observations: 119)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~labid (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.81 0.60 0.03 2.22 1.00 6487 4278
##
## ~subid_unique (Number of levels: 88)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.79 0.61 0.03 2.29 1.00 6551 3986
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 154.14 8.38 137.27 170.08 1.00 14105 5734
## CDI.z_age_months 27.43 2.93 21.44 33.06 1.00 16606 6113
## gender1 9.84 9.65 -8.35 28.43 1.00 15050 6205
## z_age_months -1.02 1.64 -5.03 1.51 1.00 9120 4618
## method1 -0.71 3.19 -9.14 4.08 1.00 6670 2125
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 91.35 6.41 79.61 104.72 1.00 15410 6275
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Calculate the marginal log likelihoods for each hypothese using Bridge sampling
marloglik_full_1.uk <- bridge_sampler(lmer.bf_full_1.uk, silent = T)
marloglik_null_1.uk <- bridge_sampler(lmer.bf_null_1.uk, silent = T)
Bayes factor
BF_full_1.uk <- bayes_factor(marloglik_full_1.uk, marloglik_null_1.uk)
BF_full_1.uk
## Estimated Bayes factor in favor of x1 over x2: 0.94293
lmer.bf_full_2.uk <- brm(vocab_nwords ~ CDI.z_age_months + gender +
z_age_months + method + IDS_pref +
IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
(1 | labid) + (1 | subid_unique),
data = data.uk,
family = gaussian,
prior = prior_uk_2,
iter = 4000,
warmup = 2000,
chains = 4,
save_pars = save_pars(all = T),
control = list(adapt_delta = .95),
seed=546)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
##
## SAMPLING FOR MODEL '8936aa99d96ef130bfa122dd46b98c39' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 4.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.41 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 1.8092 seconds (Warm-up)
## Chain 1: 1.69476 seconds (Sampling)
## Chain 1: 3.50395 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '8936aa99d96ef130bfa122dd46b98c39' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.00025 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 2.5 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 1.69324 seconds (Warm-up)
## Chain 2: 0.992679 seconds (Sampling)
## Chain 2: 2.68592 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '8936aa99d96ef130bfa122dd46b98c39' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 6.7e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.67 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 1.689 seconds (Warm-up)
## Chain 3: 1.00774 seconds (Sampling)
## Chain 3: 2.69674 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '8936aa99d96ef130bfa122dd46b98c39' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 5.8e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.58 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 1.3114 seconds (Warm-up)
## Chain 4: 0.984855 seconds (Sampling)
## Chain 4: 2.29625 seconds (Total)
## Chain 4:
summary(lmer.bf_full_2.uk)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 | labid) + (1 | subid_unique)
## Data: data.uk (Number of observations: 119)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~labid (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.81 0.61 0.04 2.26 1.00 5636 3152
##
## ~subid_unique (Number of levels: 88)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.81 0.61 0.03 2.27 1.00 6346 3454
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 154.11 8.17 137.45 169.50 1.00 8021
## CDI.z_age_months 27.41 2.93 21.50 33.06 1.00 8970
## gender1 10.00 9.79 -8.47 29.05 1.00 8842
## z_age_months -1.04 1.66 -5.25 1.43 1.00 6597
## method1 -0.69 3.10 -9.10 3.68 1.00 4440
## IDS_pref 0.55 4.80 -6.98 11.94 1.00 3006
## method1:IDS_pref 0.01 4.24 -8.69 9.16 1.00 5044
## CDI.z_age_months:IDS_pref 0.05 2.29 -5.01 5.21 1.00 5910
## z_age_months:IDS_pref -0.09 2.36 -5.13 4.70 1.00 5164
## Tail_ESS
## Intercept 5564
## CDI.z_age_months 5616
## gender1 5767
## z_age_months 4177
## method1 1560
## IDS_pref 1117
## method1:IDS_pref 1968
## CDI.z_age_months:IDS_pref 2008
## z_age_months:IDS_pref 2415
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 91.50 6.54 79.89 105.11 1.00 9123 6176
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_2.uk <- brm(vocab_nwords ~ CDI.z_age_months + gender +
z_age_months + method + IDS_pref +
IDS_pref:method + IDS_pref:CDI.z_age_months +
(1 | labid) + (1 | subid_unique),
data = data.uk,
family = gaussian,
prior = prior_uk_null_2,
iter = 4000,
warmup = 2000,
chains = 4,
save_pars = save_pars(all = T),
control = list(adapt_delta = .95),
seed =689)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
##
## SAMPLING FOR MODEL '155a5bbbeeadf50ca6dd061d19360a9f' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 5.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.55 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 1.17176 seconds (Warm-up)
## Chain 1: 0.981056 seconds (Sampling)
## Chain 1: 2.15281 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '155a5bbbeeadf50ca6dd061d19360a9f' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 5.4e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.54 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 1.21728 seconds (Warm-up)
## Chain 2: 0.949116 seconds (Sampling)
## Chain 2: 2.1664 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '155a5bbbeeadf50ca6dd061d19360a9f' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 6.1e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.61 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 1.64771 seconds (Warm-up)
## Chain 3: 0.993318 seconds (Sampling)
## Chain 3: 2.64103 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '155a5bbbeeadf50ca6dd061d19360a9f' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.000229 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 2.29 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 1.18755 seconds (Warm-up)
## Chain 4: 0.964777 seconds (Sampling)
## Chain 4: 2.15233 seconds (Total)
## Chain 4:
summary(lmer.bf_null_2.uk)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + (1 | labid) + (1 | subid_unique)
## Data: data.uk (Number of observations: 119)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~labid (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.80 0.60 0.03 2.26 1.00 7122 3819
##
## ~subid_unique (Number of levels: 88)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.80 0.61 0.03 2.26 1.00 7561 3802
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 154.24 8.33 137.43 170.45 1.00 9042
## CDI.z_age_months 27.39 2.98 21.38 33.10 1.00 10553
## gender1 10.11 9.62 -8.08 28.88 1.00 11166
## z_age_months -1.04 1.68 -5.29 1.51 1.00 6649
## method1 -0.78 3.41 -10.17 4.06 1.01 3586
## IDS_pref 0.30 4.40 -6.96 10.83 1.00 4757
## method1:IDS_pref -0.01 4.06 -8.00 8.28 1.00 4590
## CDI.z_age_months:IDS_pref 0.04 2.63 -5.12 5.61 1.01 3837
## Tail_ESS
## Intercept 5701
## CDI.z_age_months 6000
## gender1 6223
## z_age_months 3744
## method1 1188
## IDS_pref 1472
## method1:IDS_pref 1860
## CDI.z_age_months:IDS_pref 1796
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 91.34 6.45 79.76 105.22 1.00 10269 6191
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Calculate the marginal log likelihoods for each hypothese using Bridge sampling
marloglik_full_2.uk <- bridge_sampler(lmer.bf_full_2.uk, silent = T)
marloglik_null_2.uk <- bridge_sampler(lmer.bf_null_2.uk, silent = T)
Bayes factor
BF_full_2.uk <- bayes_factor(marloglik_full_2.uk, marloglik_null_2.uk)
BF_full_2.uk
## Estimated Bayes factor in favor of x1 over x2: 0.92703
lmer.bf_full_3.combined <- brm(CDI.prop ~ CDI.z_age_months + language_zone + gender +
z_age_months + method + IDS_pref + IDS_pref:language_zone +
IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
(1 + CDI.z_age_months | labid) + (1 | subid_unique),
data = data.uk_nae,
family = gaussian,
prior = prior_combined_full_info,
iter = 4000,
warmup = 2000,
chains = 4,
save_pars = save_pars(all = T),
control = list(adapt_delta = .9999),
seed=100)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
##
## SAMPLING FOR MODEL 'dc68784226044421f8a0db87fb6a9761' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000201 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.01 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 18.3068 seconds (Warm-up)
## Chain 1: 7.8792 seconds (Sampling)
## Chain 1: 26.186 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'dc68784226044421f8a0db87fb6a9761' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.000304 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 3.04 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 17.2098 seconds (Warm-up)
## Chain 2: 15.4955 seconds (Sampling)
## Chain 2: 32.7052 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'dc68784226044421f8a0db87fb6a9761' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.000152 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.52 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 16.5475 seconds (Warm-up)
## Chain 3: 7.88142 seconds (Sampling)
## Chain 3: 24.429 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'dc68784226044421f8a0db87fb6a9761' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.000151 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.51 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 16.5896 seconds (Warm-up)
## Chain 4: 15.5024 seconds (Sampling)
## Chain 4: 32.092 seconds (Total)
## Chain 4:
summary(lmer.bf_full_3.combined)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: CDI.prop ~ CDI.z_age_months + language_zone + gender + z_age_months + method + IDS_pref + IDS_pref:language_zone + IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 + CDI.z_age_months | labid) + (1 | subid_unique)
## Data: data.uk_nae (Number of observations: 424)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~labid (Number of levels: 13)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.01 0.01 0.00 0.02 1.00
## sd(CDI.z_age_months) 0.02 0.00 0.01 0.03 1.00
## cor(Intercept,CDI.z_age_months) 0.03 0.42 -0.76 0.79 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 4550 3934
## sd(CDI.z_age_months) 6399 6415
## cor(Intercept,CDI.z_age_months) 973 2167
##
## ~subid_unique (Number of levels: 302)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.01 0.01 0.00 0.03 1.00 4674 3846
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 0.32 0.01 0.29 0.34 1.00 10446
## CDI.z_age_months 0.04 0.01 0.03 0.05 1.00 5095
## language_zone1 0.08 0.02 0.05 0.11 1.00 6294
## gender1 0.03 0.01 0.01 0.05 1.00 12001
## z_age_months -0.00 0.00 -0.01 0.00 1.00 12834
## method1 -0.00 0.01 -0.02 0.02 1.00 8983
## method2 0.01 0.01 -0.02 0.04 1.00 6369
## IDS_pref 0.00 0.02 -0.03 0.04 1.00 10120
## language_zone1:IDS_pref 0.01 0.02 -0.02 0.06 1.00 5807
## method1:IDS_pref -0.00 0.02 -0.04 0.03 1.00 8979
## method2:IDS_pref -0.00 0.02 -0.04 0.04 1.00 7079
## CDI.z_age_months:IDS_pref 0.00 0.01 -0.01 0.02 1.00 14388
## z_age_months:IDS_pref -0.00 0.01 -0.02 0.01 1.00 12881
## Tail_ESS
## Intercept 6105
## CDI.z_age_months 5274
## language_zone1 5396
## gender1 5977
## z_age_months 6187
## method1 4958
## method2 3343
## IDS_pref 3866
## language_zone1:IDS_pref 2938
## method1:IDS_pref 3134
## method2:IDS_pref 3278
## CDI.z_age_months:IDS_pref 5442
## z_age_months:IDS_pref 5224
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.21 0.01 0.20 0.23 1.00 12956 5923
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_3.combined <- brm(CDI.prop ~ CDI.z_age_months + language_zone + gender +
z_age_months + method + IDS_pref +
IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
(1 + CDI.z_age_months | labid) + (1 | subid_unique),
data = data.uk_nae,
family = gaussian,
prior = prior_combined_null_info,
iter = 4000,
warmup = 2000,
chains = 4,
save_pars = save_pars(all = T),
control = list(adapt_delta = .9999),
seed =200)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
##
## SAMPLING FOR MODEL 'f762272a91b9146804973af2bbae2eb4' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000515 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.15 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 18.0606 seconds (Warm-up)
## Chain 1: 15.7231 seconds (Sampling)
## Chain 1: 33.7837 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'f762272a91b9146804973af2bbae2eb4' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.000183 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.83 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 17.0697 seconds (Warm-up)
## Chain 2: 15.3224 seconds (Sampling)
## Chain 2: 32.3921 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'f762272a91b9146804973af2bbae2eb4' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.000333 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 3.33 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 18.7739 seconds (Warm-up)
## Chain 3: 8.02287 seconds (Sampling)
## Chain 3: 26.7968 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'f762272a91b9146804973af2bbae2eb4' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.000525 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 5.25 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 18.542 seconds (Warm-up)
## Chain 4: 7.94283 seconds (Sampling)
## Chain 4: 26.4848 seconds (Total)
## Chain 4:
summary(lmer.bf_null_3.combined)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: CDI.prop ~ CDI.z_age_months + language_zone + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 + CDI.z_age_months | labid) + (1 | subid_unique)
## Data: data.uk_nae (Number of observations: 424)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~labid (Number of levels: 13)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.01 0.01 0.00 0.02 1.00
## sd(CDI.z_age_months) 0.02 0.00 0.01 0.03 1.00
## cor(Intercept,CDI.z_age_months) 0.01 0.43 -0.79 0.79 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 5150 4137
## sd(CDI.z_age_months) 7566 6438
## cor(Intercept,CDI.z_age_months) 1200 2826
##
## ~subid_unique (Number of levels: 302)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.01 0.01 0.00 0.03 1.00 4593 3446
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 0.31 0.01 0.29 0.34 1.00 10940
## CDI.z_age_months 0.04 0.01 0.03 0.05 1.00 6321
## language_zone1 0.08 0.02 0.05 0.11 1.00 7368
## gender1 0.03 0.01 0.01 0.05 1.00 15525
## z_age_months -0.00 0.00 -0.01 0.00 1.00 14549
## method1 -0.00 0.01 -0.02 0.02 1.00 11072
## method2 0.01 0.02 -0.02 0.04 1.00 7455
## IDS_pref 0.00 0.02 -0.03 0.03 1.00 11181
## method1:IDS_pref -0.00 0.02 -0.04 0.03 1.00 9539
## method2:IDS_pref 0.00 0.02 -0.03 0.04 1.00 9235
## CDI.z_age_months:IDS_pref 0.00 0.01 -0.01 0.02 1.00 15202
## z_age_months:IDS_pref -0.00 0.01 -0.02 0.01 1.00 12744
## Tail_ESS
## Intercept 6379
## CDI.z_age_months 5664
## language_zone1 5211
## gender1 5569
## z_age_months 6425
## method1 5420
## method2 3905
## IDS_pref 4151
## method1:IDS_pref 3629
## method2:IDS_pref 3787
## CDI.z_age_months:IDS_pref 4821
## z_age_months:IDS_pref 5179
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.21 0.01 0.20 0.23 1.00 13563 6219
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Calculate the marginal log likelihoods for each hypothese using Bridge sampling
marloglik_full_3.combined <- bridge_sampler(lmer.bf_full_3.combined, silent = T)
marloglik_null_3.combined <- bridge_sampler(lmer.bf_null_3.combined, silent = T)
Bayes factor
BF_full_3.combined <- bayes_factor(marloglik_full_3.combined, marloglik_null_3.combined)
BF_full_3.combined
## Estimated Bayes factor in favor of x1 over x2: 0.06768
```